{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solvers for Linear Algebraic Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have implemented multigrid solvers for linear algebraic systems arising from various finite element methods. Here we briefly present the usage for a symmetric and positive definite matrix equation `Ax = b`.\n",
    "\n",
    "    x = A\\b;\n",
    "    x = amg(A,b);\n",
    "    x = mg(A,b,elem);\n",
    "    \n",
    "Backslash `\\` is the build-in direct solver of MATLAB. It is suitable for small size problems. `x = amg(A,b)` implements algebraic multigrid solver. To acheive multigrid efficiency, a hierarchical 'grids' is generated from the graph of `A`. When the mesh is avaiable, `x = mg(A,b,elem)` implements geometric multigrid solvers. Inside `mg`, an coarsening algorithm is applied to the mesh given by `elem`. \n",
    "\n",
    "More options is allowed in `mg` and `amg`. Try `help mg` and `help amg` for possible options including tolerance, V or W-cycles, number of smoothings steps, and print level etc.\n",
    "\n",
    "Here we include several examples for discrete Poisson matrices. Solvers for other finite element methods and other equations can be found\n",
    "- [List of Examples](solverlist.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: 2-D Linear Element "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Direct solver\n",
      "Elapsed time is 1.293569 seconds.\n",
      "\n",
      " Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:   263169,  #nnz:  1303561, smoothing: (1,1), iter: 10,   err = 1.35e-09,   time =  1.1 s\n",
      "Elapsed time is 1.186739 seconds.\n",
      "\n",
      " Algebraic Multigrid W-cycle Preconditioner with Conjugate Gradient Method\n",
      "  nnz/N: 4.96,   level:  6,   coarse grid 169,   nnz/Nc 9.38\n",
      "#dof:  263169,    iter: 14,   err = 5.5672e-09,   time = 4.65 s\n",
      " \n",
      "Elapsed time is 4.334722 seconds.\n",
      "Difference between direct and mg, amg solvers 1.4e-09, 7.8e-08 \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4gsWBCE4Nm4s2gAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1Ob3YtMjAxOCAyMDozMzo1Ni9aGsoAAA1Z\nSURBVHic7do9jh3HGYXhvoZyewEkwFyJUjKZCRUqMLQOAbMBbUCwl+DUjBQqnEnElAlzAjQUE9zA\ndTAEf4Z37m9311d1nmcFJ6oX3VWb7XY7AUBrf2s9AACmSZAAKEKQAChBkAAoQZAAKEGQAChBkAAo\nQZAAKEGQAChBkAAoQZAAKEGQAChBkAAoQZAAKEGQAChBkAAoQZAAKEGQAChBkAAoQZAAKEGQAChB\nkAAoQZAAKEGQAChBkAAoQZAAKEGQAChBkAAoQZAAKOG71gPGd3t7e3d313oF0MzV1dX19XXrFR3Y\nbLfb1hsGt9lsrm+et15BG+/ffXj98s2z50+fvXjSegtt3P72apomJ+0xBGlxm83m179uWq+ggffv\nPvznn//9x5O/P3vx5PrmRes5NPD7L39M0/T65Rsn7THcIcEi7mv0079+9G0U675GP/37x9ZDuiFI\nML8vavS09RbaUKMzCBLMTI1Qo/MIEsxJjVCjswkSzEaNUKNLCBLMQ41QowsJEsxAjVCjywkSXEqN\nUKNZCBJcRI1Qo7kIEpxPjVCjGQkSnEmNUKN5CRKcQ41Qo9kJEpxMjVCjJQgSnEaNUKOFCBKcQI1Q\no+UIEhxLjVCjRQkSHEWNUKOlCRIcpkao0QoECQ5QI9RoHYIE+6gRarQaQYJHqRFqtCZBgt3UCDVa\nmSDBDmqEGq1PkOAhNUKNmhAk+IoaoUatCBJ8pkaoUUOCBB+pEWrUliDBNKkRalSAIIEaoUYlCBLp\n1Ag1KkKQiKZGqFEdgkQuNUKNShEkQqkRalSNIJFIjVCjggSJOGqEGtUkSGRRI9SoLEEiiBqhRpUJ\nEinUCDUqTpCIoEaoUX2CxPjUCDXqgiAxODVCjXohSIxMjVCjjggSw1Ij1KgvgsSY1Ag16o4gMSA1\nQo16JEiMRo1Qo04JEkNRI9SoX4LEONQINeqaIDEINUKNeidIjECNUKMBCBLdUyPUaAyCRN/UCDUa\nhiDRMTVCjUYiSPRKjVCjwQgSXVIj1Gg8gkR/1Ag1GpIg0Rk1Qo1GJUj0RI1Qo4EJEt1QI9RobIJE\nH9QINRqeINEBNUKNEggS1akRahRCkChNjVCjHIJEXWqEGkURJIpSI9QojSBRkRqhRoEEiXLUCDXK\nJEjUokaoUSxBohA1Qo2SCRJVqBFqFE6QKEGNUCMEifbUCDViEiSaUyPUiHuCREtqhBrxiSDRjBqh\nRnxJkGhDjVAjHhAkGlAj1IhvCRJrUyPUiJ0EiVWpEWrEYwSJ9agRasQegsRK1Ag1Yj9BYg1qhBpx\nkCCxODVCjTiGILEsNUKNOJIgsSA1Qo04niCxFDVCjTiJILEINUKNOJUgMT81Qo04gyAxMzVCjTiP\nIDEnNUKNOJsgMRs1Qo24hCAxDzVCjbiQIDEDNUKNuJwgcSk1Qo2YhSBxETVCjZiLIHE+NUKNmJEg\ncSY1Qo2YlyBxDjVCjZidIHEyNUKNWIIgcRo1Qo1YiCBxAjVCjViOIHEsNUKNWJQgcRQ1Qo1YmiBx\nmBqhRqxAkDhAjVAj1iFI7KNGqBGrESQepUaoEWsSJHZTI9SIlQkSO6gRasT6BImH1Ag1oglB4itq\nhBrRiiDxmRqhRjQkSHykRqgRbQkS06RGqBEFCBJqhBpRgiClUyPUiCIEKZoaoUbUIUi51Ag1ohRB\nCqVGqBHVCFIiNUKNKEiQ4qgRakRNgpRFjVAjyhKkIGqEGlGZIKVQI9SI4gQpghqhRtQnSONTI9SI\nLgjS4NQINaIXgjQyNUKN6IggDUuNUCP6IkhjUiPUiO4I0oDUCDWiR4I0GjVCjeiUIA1FjVAj+iVI\n41Aj1IiuCdIg1Ag1oneCNAI1Qo0YgCB1T41QI8YgSH1TI9SIYQhSx9QINWIkgtQrNUKNGIwgdUmN\nUCPGI0j9USPUiCEJUmfUCDViVILUEzVCjRiYIHVDjVAjxiZIfVAj1IjhCVIH1Ag1IoEgVadGqBEh\nBKk0NUKNyCFIdakRakQUQSpKjVAj0ghSRWqEGhFIkMpRI9SITIJUixqhRsQSpELUCDUimSBVoUao\nEeEEqQQ1Qo1AkNpTI9QIJkFqTo1QI7gnSC2pEWoEnwhSM2qEGsGXBKkNNUKN4AFBakCNUCP4liCt\nTY1QI9hJkFalRqgRPEaQ1qNGqBHsIUgrUSPUCPYTpDWoEWoEBwnSGtQo3OuXbyY1gkM22+229YbB\nbTab1hOAxpy0x/iu9YAIv/5103oCbfz+yx9vX7374efvr29etN5CA/e/69+/+9B6SB/8soOl3N8b\n/fDz962H0Many+PWQ7ohSLAIrxjCecp0BkGC+alRODU6jyDBzNQonBqdTZBgTmoUTo0uIUgwGzUK\np0YXEiSYhxqFU6PLCRLMQI3CqdEsBAkupUbh1GguggQXUaNwajQjQYLzqVE4NZqXIMGZ1CicGs1O\nkOAcahROjZYgSHAyNQqnRgsRJDiNGoVTo+UIEpxAjcKp0aIECY6lRuHUaGmCBEdRo3BqtAJBgsPU\nKJwarUOQ4AA1CqdGqxEk2EeNwqnRmgQJHqVG4dRoZYIEu6lRODVanyDBDmoUTo2aECR4SI3CqVEr\nggRfUaNwatSQIMFnahROjdoSJPhIjcKpUXOCBNOkRvHUqAJBAjVKp0ZFCBLp1CicGtUhSERTo3Bq\nVIogkUuNwqlRNYJEKDUKp0YFCRKJ1CicGtUkSMRRo3BqVJYgkUWNwqlRZYJEEDUKp0bFCRIp1Cic\nGtUnSERQo3Bq1AVBYnxqFE6NeiFIDE6NwqlRRwSJkalRODXqiyAxLDUKp0bdESTGpEbh1KhHgsSA\n1CicGnVKkBiNGoVTo34JEkNRo3Bq1DVBYhxqFE6NeidIDEKNwqnRAASJEahRODUagyDRPTUKp0bD\nECT6pkbh1GgkgkTH1CicGg1GkOiVGoVTo/EIEl1So3BqNCRBoj9qFE6NRiVIdEaNwqnRwASJnqhR\nODUamyDRDTUKp0bDEyT6oEbh1CiBINEBNQqnRiEEierUKJwa5RAkSlOjcGoURZCoS43CqVEaQaIo\nNQqnRoEEiYrUKJwaZRIkylGjcGoUS5CoRY3CqVEyQaIQNQqnRuEEiSrUKJwaIUiUoEbh1IhJkKhA\njcKpEfcEicbUKJwa8Ykg0ZIahVMjviRINKNG4dSIBwSJNtQonBrxLUGiATUKp0bsJEisTY3CqRGP\nESRWpUbh1Ig9BIn1qFE4NWI/QWIlahROjThIkFiDGoVTI44hSCxOjcKpEUcSJJalRuHUiOMJEgtS\no3BqxEkEiaWoUTg14lSCxCLUKJwacQZBYn5qFE6NOI8gMTM1CqdGnE2QmJMahVMjLiFIzEaNwqkR\nFxIk5qFG4dSIywkSM1CjcGrELASJS6lRODViLoLERdQonBoxI0HifGoUTo2YlyBxJjUKp0bMTpA4\nhxqFUyOWIEicTI3CqRELESROo0bh1IjlCBInUKNwasSiBIljqVE4NWJpgsRR1CicGrECQeIwNQqn\nRqxDkDhAjcKpEasRJPZRo3BqxJoEiUepUTg1YmWCxG5qFE6NWJ8gsYMahVMjmhAkHlKjcGpEK4LE\nV9QonBrRkCDxmRqFUyPaEiQ+UqNwakRzgsQ0qVE8NaICQUKN0qkRRQhSOjUKp0bUIUjR1CicGlGK\nIOVSo3BqRDWCFEqNwqkRBQlSIjUKp0bUJEhx1CicGlGWIGVRo3BqRGWCFESNwqkRxQlSCjUKp0bU\nJ0gR1CicGtEFQRqfGoVTI3ohSINTo3BqREcEaWRqFE6N6IsgDUuNwqkR3RGkMalRODWiR4I0IDUK\np0Z0SpBGo0bh1Ih+CdJQ1CicGtE1QRqHGoVTI3onSINQo3BqxAAEaQRqFE6NGIMgdU+NwqkRwxCk\nvqlRODViJILUMTUKp0YMRpB6pUbh1IjxCFKX1CicGjEkQeqPGoVTI0YlSJ1Ro3BqxMAEqSdqFE6N\nGJsgdUONwqkRwxOkPqhRODUigSB1QI3CqREhBKk6NQqnRuQQpNLUKJwaEUWQ6lKjcGpEGkEqSo3C\nqRGBBKkiNQqnRmQSpHLUKJwaEUuQalGjcGpEMkEqRI3CqRHhBKkKNQqnRiBIJahRODWCSZAqUKNw\nagT3BKkxNQqnRvCJILWkRuHUCL4kSM2oUTg1ggcEqQ01CqdG8C1BakCNwqkR7CRIa1OjcGoEjxGk\nValRODWCPQRpPWoUTo1gP0FaiRqFUyM4SJDWoEbh1AiOsdlut603DG6z2UzTdH3zvPUQ2nj75//e\nvnr37PnTZy+etN5CG7e/vXLSHkOQFnd7e3t3d9d6BdDM1dXV9fV16xUdECQASnCHBEAJggRACYIE\nQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRA\nCYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJ\nggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACYIEQAmCBEAJggRACf8H1Hg8\nqOD7E6gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "% setting\n",
    "mesh.shape = 'square';\n",
    "mesh.type = 'uniform';\n",
    "mesh.size = 2e5;\n",
    "pde = 'Poisson';\n",
    "fem = 'P1';\n",
    "% get the matrix\n",
    "[eqn,T] = getfemmatrix(mesh,pde,fem);\n",
    "% compare solvers\n",
    "tic; disp('Direct solver'); x1 = eqn.A\\eqn.b; toc;\n",
    "tic; x2 = mg(eqn.A,eqn.b,T.elem); toc;\n",
    "tic; x3 = amg(eqn.A,eqn.b); toc;\n",
    "format shorte\n",
    "fprintf('Difference between direct and mg, amg solvers %0.2g, %0.2g \\n',...\n",
    "         norm(x1-x2)/norm(eqn.b),norm(x1-x3)/norm(eqn.b));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For problem size of $2.6 \\times 10^5$, `mg` ties with direct solver `\\`. But `amg` is aroud 3-4 times slover."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: 2-D Adaptive Meshes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Direct solver\n",
      "Elapsed time is 3.347438 seconds.\n",
      "\n",
      " Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:   738561,  #nnz:  3682043, smoothing: (1,1), iter: 10,   err = 6.29e-09,   time =  3.1 s\n",
      "Elapsed time is 2.843210 seconds.\n",
      "\n",
      " Algebraic Multigrid W-cycle Preconditioner with Conjugate Gradient Method\n",
      "  nnz/N: 4.99,   level:  6,   coarse grid 453,   nnz/Nc 9.72\n",
      "#dof:  738561,    iter: 15,   err = 2.9982e-09,   time = 14.9 s\n",
      " \n",
      "Elapsed time is 13.711571 seconds.\n",
      "Difference between direct and mg, amg solvers 1.2e-08, 6e-08 \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4gsWBC4GcJctvgAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1Ob3YtMjAxOCAyMDo0NjowNtfTwq0AACAA\nSURBVHic7d2/bmNHloDxq0XHFOOFLsBcGygzqGSUTYcN2KtOG1h0Lq8ewH4BwcqddNqEB1CobOiE\nRGcMhjmBK2xM6QW0wbFpDUWR99atP6dOfb9oEjXLQg+/Jm9VnaPn5+cKAIDU/iP1AgAAqCqCBABQ\ngiABAFQgSAAAFQgSAEAFggQAUIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQ\ngSABAFQgSAAAFQgSAEAFggQAUIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQ\ngSABAFQgSAAAFQgSAECFd6kXYN90Ov39999TrwJAMn/7298uLi5SryIDR8/Pz6nXYNzR0dGwHpxd\nnqZeSBrr5mkxWRb7G5D//Kqqzi5Ph/Ug9XISWM0eVvOm2P/8qqqmN/OqqninbYNPSDHIe/HF9Xnq\nhSSwmjWreTMa11V5v4F18/Tlh6+jcT2sB+vm6eL6vLQ35dWsWc0eRuP67PJ0dF6nXk4Cd1f3Z5en\n8o8SHMQzpBjkjXh6M0u9kDSGJ8cX1+fr5qmo34DU6MMv70fnJ8N68OH2/d3V/bp5Sr2ueFazZnoz\n//SPy9QLSebu6r6qqg+371MvJBsEKZLSm1QPimrSixr98bGgtCZRI2rkgCDFQ5MKadLrGolymkSN\nqJEbghQVTTLfpLdqJEpoEjWiRs4IUmw0yXCT9tdI2G4SNaJGfRCkBGiSySa1qZGw2iRqRI16Ikhp\n0CRjTWpfI2GvSdSIGvVHkJKhSWaa1LVGwlKTqBE18oIgpUSTDDTJrUbCRpOoETXyhSAlRpOyblKf\nGoncm0SNqJFHBCk9mpRpk/rXSOTbJGpEjfwiSCrQpOya5KtGIscmUSNq5B1B0oImZdQkvzUSeTWJ\nGlGjEAiSIjQpiyaFqJHIpUnUiBoFQpB0oUnKmxSuRkJ/k6gRNQqHIKlDk9Q2KXSNhOYmUSNqFBRB\n0ogmKWxSnBoJnU2iRtQoNIKkFE1S1aSYNRLamkSNqFEEBEkvmqSkSfFrJPQ0iRpRozgIkmo0KXmT\nUtVIaGgSNaJG0RAk7WhSwialrZFI2yRqRI1iIkgZoElJmqShRiJVk6gRNYqMIOWBJkVukp4aifhN\nokbUKD6ClA2aFK1J2mokYjaJGlGjJAhSTmhShCbprJGI0yRqRI1SIUiZoUlBm6S5RiJ0k6gRNUqI\nIOWHJgVqkv4aiXBNokbUKC2ClCWa5L1JudRIhGgSNaJGyRGkXNEkj03Kq0bCb5OoETXSgCBljCZ5\naVKONRK+mkSNqJESBClvNKlnk/KtkejfJGpEjfQgSNmjSc5Nyr1Gok+TqBE1UoUgWUCTHJpko0bC\nrUnUiBppQ5CMoEmdmmSpRqJrk6gRNVKIINlBk1o2yV6NRPsmUSNqpBNBMoUmHWyS1RqJNk2iRtRI\nLYJkDU3a0yTbNRL7m0SNqJFmBMkgmrSzSSXUSLzVJGpEjZQjSDbRpK0mlVMj8bpJ1Iga6UeQzKJJ\nmyaVViPxsknUiBplgSBZRpOkSbff/VpajcSmSdSookY5IEjGFd6kqqoWk+WwHqzmTeqFpLFuHlfz\nZv3wGG32uTbUKCMEyb5imyTf1J1dng5PjuPMPtdGvqk7uzw9uzyNM/tcG2qUF4JUhAKb9PK5UZzZ\n59q8fG40rI8jzD7XhhplhyCVoqgmvd7FUFqTXu9iCD37XBtqlCOCVJBCmvTWnrpymvTWnrpymkSN\nMkWQymK+Sft3eJfQpP07vEtoEjXKF0EqjuEmtTlvZLtJbc4b2W4SNcoaQSqRySa1P/1qtUntT79a\nbRI1yh1BKpSxJnW9i8Fek7rexWCvSdTIAIJULjNNcrsZyFKT3G4GstQkamQDQSqagSb1uafORpP6\n3FNno0nUyAyCVLqsm9T/1tTcm9T/1tTcm0SNLCFIyLVJvu7wzrdJvu7wzrdJ1MgYgoSqyrBJfidK\n5NgkvxMlcmwSNbKHIOEPGTUpxHyjvJoUYr5RXk2iRiYRJPwliyaFm7aXS5PCTdvLpUnUyCqChH+j\nvEmhZ7/qb1Lo2a/6m0SNDCNI2Ka2SXEmkWtuUpxJ5JqbRI1sI0jYQWGT4tRI6GxSnBoJnU2iRuYR\nJOymqkkxayS0NSlmjYS2JlGjEhAkvElJk+LXSOhpUvwaCT1NokaFIEjYJ3mTUtVIaGhSqhoJDU2i\nRuUgSDggYZPS1kikbVLaGom0TaJGRSFIOCxJkzTUSKRqkoYaiVRNokalIUhoJXKT9NRIxG+SnhqJ\n+E2iRgUiSGgrWpO01UjEbJK2GomYTaJGZSJI6CBCk3TWSMRpks4aiThNokbFIkjoJmiTNNdIhG6S\n5hqJ0E2iRiUjSOgsUJP010iEa5L+GolwTaJGhSNIcOG9SbnUSIRoUi41EiGaRI1AkODIY5PyqpHw\n26S8aiT8NokaoSJI6MNLk3KskfDVpBxrJHw1iRpBECT00rNJ+dZI9G9SvjUS/ZtEjbBBkNCXc5Ny\nr5Ho06TcayT6NIka4SWCBA8cmmSjRsKtSTZqJNyaRI2whSDBj05NslQj0bVJlmokujaJGuE1ggRv\nWjbJXo1E+ybZq5Fo3yRqhJ0IEnw62CSrNRJtmmS1RqJNk6gR3kKQ4NmeJtmukdjfJNs1EvubRI2w\nB0GCfzubVEKNxFtNKqFG4q0mUSPsR5AQxFaTyqmReN2kcmokXjeJGuEggoRQ/mzSvKqqomokNk1a\nTJar2UNRNRJ/NenhUf4aUCPsd/T8/Jx6DcYdHR2NxgW9EW9ZzZuqqob1YHhynGQB64fHqqoSvrp8\nREj4d4DfQHKrecM7bRvvUi+gCKt58+m3sv51LNYPT/J+NDw5vrgeJ1nDat4sJsuqqpIsYDVv5MNB\nqgVUVTW9mY/OT1L1YHozlyCdfTwdngySrCGtux/vUy8hGwQphg+37+9+vL/69jn1QqJaN093P96f\nXZ6uZg+j85PVvJEv8WIv4+FpeHKcZAGrWbOaPVxcj+UdOdVvYDhZDuvjJN+X3l3dD+vBqKrPPp4u\nvi4/3L4f1mU16cv3kw+/vP/ywyT1QvLAM6QYzi5PL67Pb7/7NfVC4vlrF8O4rqLMPt9jWA/OLv8r\n8gJe7mKIM/tcm5e7GIYnMWafa/Pl+8nF9bioR6c9EaRIimrSzj11RTXp9Z660pr0ek9d6Nnn2lAj\nBwQpnkKatGeHdyFNemuHdzlNemuHdzlNokZuCFJU5pt08LyR+SbtP29UQpP2nzcqoUnUyBlBis1w\nk1qefjXcpDanX203qc3pV9tNokZ9EKQETDap010MJpvU/i4Gq01qfxeD1SZRo54IUhrGmuRwM5Cx\nJnW9Gchek7reDGSvSdSoP4KUjJkmOd9TZ6ZJbvfUWWqS2z11lppEjbwgSCkZaFLPW1MNNKnPrak2\nmtTn1lQbTaJGvhCkxLJukpc7vLNuUv87vHNvUv87vHNvEjXyiCCll2mTPE6UyLRJviZK5NskXxMl\n8m0SNfKLIKmQXZO8zzfKrkl+5xvl2CS/841ybBI18o4gaZFRkwJN28uoSSGm7eXVpBDT9vJqEjUK\ngSApkkWTgs5+zaJJ4Wa/5tKkcLNfc2kSNQqEIOmivEkRJpErb1LoSeT6mxR6Ern+JlGjcAiSOmqb\nFKFGQm2TQtdoswC1TQpdI6G5SdQoKIKkkcImRauRUNikODXaLEBhk+LUSOhsEjUKjSAppapJkWsk\nVDUpZo02C1DVpJg1EtqaRI0iIEh6KWlSkhoJJU26u7qPXKPNApQ0KX6NhJ4mUaM4CJJqyZuUsEYi\neZOG9fFishydn6RaQPImpaqR0NAkahQNQdIuYZOS10gkbNJq1iy+Lof1INUCqtRNSlsjkbZJ1Cgm\ngpSBJE1SUiORpEkvnxvFmX3+llRN0lAjkapJ1CgygpSHyE1SVSMRuUlbuxgizD7fL36T9NRIxG8S\nNYqPIGUjWpMU1khEa9LOPXVFNUlbjUTMJlGjJAhSTiI0SW2NRIQm7dnhXUiTdNZIxGkSNUqFIGUm\naJOU10gEbdLB80bmm6S5RiJ0k6hRQgQpP4GalEWNRKAmtTz9arhJ+mskwjWJGqVFkLLkvUkZ1Uh4\nb1KnuxhMNimXGokQTaJGyRGkXHlsUnY1Eh6b5HAzkLEm5VUj4bdJ1EgDgpQxL03KtEbCS5Oc76kz\n06QcayR8NYkaKUGQ8tazSVnXSPRsUs9bUw00Kd8aif5NokZ6EKTsOTfJQI2Ec5O83OGddZNyr5Ho\n0yRqpApBssChSWZqJBya5HGiRKZNslEj4dYkaqQNQTKiU5OM1Uh0apL3+UbZNclSjUTXJlEjhQiS\nHS2bZLJGomWTAk3by6hJ9mok2jeJGulEkEw52CTDNRIHmxR09msWTbJaI9GmSdRILYJkzZ4mma+R\n2NOkCJPIlTfJdo3E/iZRI80IkkE7m1RIjcTOJkWokVDbpBJqJN5qEjVSjiDZtNWkomoktpoUrUZC\nYZPKqZF43SRqpB9BMmvTpAJrJKRJi6/L1TxqjcSmSYvJMubrvlyANGk1byRL5dRIvGwSNcrC0fPz\nc+o1GHd0dPTz/12nevXpzWx6Mx/Wg1QLWDdPCV9dFlBVVZ819PlP6P/qPW0+HxT7d0B+Ax9u359d\nnqZaw8//ecM7bRvvUi8AAa2bp8VkeXZ5upgsr759TrSGxy8/TEbjOsk/z1fz5u7qflgPRuNaPjA5\n+PLD13Xz9Om3y2F93PVnF5N/TW/mVVX1WUAf05uZfERLtYCqqu6u7s8+no7GaT6dyKsvvi5H4zrt\nv41wEEEy6+U3daPz+ssPX5M0ad08jsb16PxkejOL3KTVrFl8XX64fb+aNcN6sJj8y/kd+erb5y8/\nfP3028dO72iLyXI1e7i4HstHhD4LcCPPjc4uT4f1YN08xV/AxvBkkCQGm2/qRuP67ur+w+17mqQZ\nz5Bs2npuFGH2+X4X1+fDeiDvj3Fs7WLoeQfrsB58+u2jfFRq+SOLyXLxdSkLSLLH4eUuhmF9HGH2\nuTYvnxvFmX2OngiSQTt3MRTVpJ176mI26WWNNj8es0mv99SFnn2uzetdDDRJP4JkzZ49dYU0ac8O\n7zhNel2jzY/HadJbO7zLadJbe+poknIEyZSDO7zNN+ngeaPQTXqrRpsfD92k/eeNSmjS/h3eNEkz\ngmRHy/NGhpvU8vRruCbtr9Hmx8M1qc3pV9tNanPeiCapRZCM6HT61WSTOt3FEKJJbWq0+fEQTWp/\nF4PVJrU//UqTdCJIFjjcxWCsSQ43A/ltUvsabX7cb5O63gxkr0ld72KgSQoRpOw53wxkpknO99T5\natL0ZtapRpsf99Ukt3vqLDXJ7WYgmqQNQcpbz3vqDDSp562p/Zt0dnk6vZm7nfn10qQ+t6baaFKf\ne+pokioEKWNebk3Nukle7vDu0yS5i0HucXB7R+vZpP53eOfepP63ptIkPQhSrjze4Z1pkzxOlHBr\n0ua5kcM9Di85N8nXRIl8m+TrDm+apARBypL3iRLZNcn7fKOuTdraxRC/SX7nG+XYJL8TJWiSBgQp\nP4HmG2XUpEDT9to3aeeeuphNCjFtL68mhZhvRJOSI0iZCTptL4smBZ392qZJe3Z4x2lSuNmvuTQp\n3LQ9mpQWQcpJhNmvypsUYRL5/iYdPG8UukmhJ5Hrb1Lo2a80KSGClI1ok8jVNilCjTYLqHYloeXp\n13BNCl2jzQLUNinOJHKalApBykO0GgmFTYpWo80Cqn9PQqe7GEI0KU6NNgtQ2KQ4NRI0KQmClIHI\nNRKqmhS5RpsFVH8moevNQJXvJsWs0WYBqpoUs0aCJsXHCHPtktRInF2eVlV1+92vSWafV1V1cX0+\nvZl9+X6yfni8uD5fTJYOf8hq1qybJ7efHdbH05uZvCm7LUDucTi7PHVbQPXn2Vt5c3T7E5xJk6Y3\nM/nPj/zqL8Wvkdg0idnncRAk1RLWSCRv0mhcy1v5unl0+xPWzdNq3gzrQc83FLcFbK5erarKbQHr\n5mndPF2cjx1+tj8NTUpVI0GTYiJIeiWvkUjYJPmm7urbZ/mY4vYRYVFLz57OLk+7/iYXk+Xw5Lj6\n85fQ9R153TzdfvfraFwP68Fq3ny6/tj1He3u6n7z6SpVEtI2KW2NBE2KhmdISimpkUjyPOnlc6Oe\nd7DKG8r0Zr6aNe1/6uVzI4e7haRGn367HJ2fuD1P2jw3ijb7/C2pnidpqJHgeVIcBEkjVTUSkZv0\nehdD5Ca93sXQqUkvalRvFtCpSVu7GApskp4aCZoUAUFSR2GNRLQmvbWnLlqT3tpT17JJr2u0WUDL\nJu3cU1dUk7TVSNCk0AiSLmprJCI0af8O7whN2r/D+2CT3qrRZgEHm7Rnh3chTdJZI0GTgiJIiiiv\nkQjapDbnjYI2qc15oz1N2l+jzQL2NOngeSPzTdJcI0GTwiFIWmRRIxGoSe1PvwZqUvvTrzub1KZG\nmwXsbFLL06+Gm6S/RoImBUKQVMioRsJ7k7rexeC9SV3vYthqUvsabRaw1aROdzGYbFIuNRI0KQSC\nlF52NRIem+R2M5DHJjncDFS9aFLXGm0WsGmSw81AxpqUV40ETfKOICWWaY2Elyb1uafOV5McarRZ\ngFuNNgv49NtH+QU6HPs106QcayRokl8EKaWsayR6Nqn/rak9m7SaN5uVOPy43JJ3dnm6+XO6mt7M\n5MejzT73q3+T8q2RoEkeEaRkDNRIODfJ1x3ezk3afFPncI9D9eK5kXy4cXhH3nxTF232eQh9mpR7\njQRN8oUgpWGmRsKhSX4nSjg06eVzI4e7hba+qXO4W+jlc6M4s8/DcWuSjRoJmuQFQUrAWI1EpyaF\nmG/UqUmvdzF0atLO50admvR6F0NpTbJUI0GT+iNIsZmskWjZpHDT9lo26a09dS2btGcXQ8smvbWn\nrpwm2auRoEk9EaSoDNdIHGxS6NmvB5u0f4f3wSYd3FN3sEn7d3iX0CSrNRI0qQ+CFI/5Gok9TYoz\niXxPk9qcN9rTpJY7vPc0qc15I9tNsl0jQZOcEaRICqmR2NmkODUSO5vU/vTrziZ1Om+0s0ntT79a\nbVIJNRI0yQ1BiqGoGomtJsWskdhqUte7GLaa5HD6datJXe9isNekcmokaJKDo+fn59RrMO7o6GhY\nD4qq0cZisry7uh/Wg+HJccwabUxvZjL/220BcqmPHHp1u4ths4DRuHa4i0H+KVNVldv8+HXztJj8\nazFZXlyfyxT2yCRIi8my2P8LbP4K8U7bBkEK7ujoKPUSAKR0cXHxz3/+M/UqMvAu9QKKIP86HtaD\n1AuJbTVr7n68XzdPw3rg9m/8/haT5fRm5vYBRdx+9+vw5Njt6yb51/HV//zvTz/95PbqQDl4hhRD\nmV8ly3OjD7+8H43rOLPP3zIa133uu6uqyu1uIXnFJN+VATkiSDEU+HhzaxdDhNnn+0WYfb7FYaIE\nUDiCFElRTdq5p66oJlEjwAFBiqeQJu3Z4V1Ik6gR4IYgRWW+SQfPG5lvEjUCnBGk2Aw3qeXpV8NN\nokZAHwQpAZNN6nQXg8kmUSOgJ4KUhrEmOdwMZKxJ1AjojyAlY6ZJzvfUmWnSl+8nFTUCeiNIKRlo\nUs9bUw00aVgPVvOG069AfwQpsayb5OUO76ybJD919e2zwz0OALYQpPQybZLHiRKZNmnz3MjhHgcA\nrxEkFbJrkvf5Rtk1aWsXA00C+iNIWmTUpEDT9jJq0s49dTQJ6IkgKZJFk4LOfs2iSXt2eNMkoA+C\npIvyJkWYRK68SQfPG9EkwBlBUkdtkyLUSKhtUsvTrzQJcEOQNFLYpGg1Egqb1OkuBpoEOCBISqlq\nUuQaCVVNcrgZ6EWTHkItEbDlXeoF4E2bJslJl1TLSFIjIdcf3H7369W3z/Ffvaqqi+vz6c1sPXla\nNMuL6/H0Ztb1TxjWg8VkWf09xOoAawiSasmblLBGInmT1s3TaFyvHx6H9bHDj7v9FFAmgqRdwiYl\nr5FI2CT5pu7TPy6nN7PVrHG7PnXdPPpeF2ATz5AykOR5kpIaiSTPk14+N+p5ByuANghSHiI3SVWN\nROQmvd7FQJOA0AhSNqI1SWGNRLQmvbWnjiYBQRGknERoktoaiQhN2r/DmyYB4RCkzARtkvIaiaBN\nanPeiCYBgRCk/ARqUhY1EoGa1P70K00CQiBIWfLepIxqJLw3qetdDDQJ8I4g5cpjk7KrkfDYJIeb\ngSqaBPhGkDLmpUmZ1kh4aZJbjQRNAjwiSHnr2aSsayR6NqlPjQRNAnwhSNlzbpKBGgnnJvWvkaBJ\ngBcEyQKHJpmpkXBokq8aCZoE9EeQjOjUJGM1Ep2a5LdGgiYBPREkO1o2yWSNRMsmhaiRoElAHwTJ\nlINNMlwjcbBJ4WokaBLgjCBZs6dJ5msk9jQpdI0ETQLcECSDdjapkBqJnU2KUyNBkwAHBMmmrSYV\nVSOx1aSYNRI0CeiKEeZmbZp09vF0ejP79NvHmANnN9YPT1VVJXnp0bg+uzyVJIzG9Yfb95GXcXb5\nX9Ob2fRmfvHz32O+LpApgmTZsB6cfTy9u7of1oMvP3xNtYx185Tw1f9Yw8Nj8jUA2I8gWbaaNYuv\ny6tvn++u7j/cvh/WgyRrSPht4d3V/apq5H9fffucZA3Tm1mS1wWywzMkszYliDb7XBv5su7i+nw0\nruPMPgfQB0GyaetzSYFN2trFEGH2OYCeCJJBO78lK6pJO/fU0SRAOYJkzZ5nNoU0ac8Ob5oEaEaQ\nTDm4g8B8kw6eN6JJgFoEyY6W+9kMN6nl6VeaBOhEkIzotLvaZJM63cVAkwCFCJIFDmd9jDXJ4WYg\nmgRoQ5Cy53zy1EyTnO+po0mAKgQpbz3vQTDQpJ63ptIkQA+ClDEvt/Jk3SQvd3jTJEAJgpQrj3fE\nZdokjxMlaBKgAUHKkvcbS7Nrkvf5RjQJSI4g5SfQ/dkZNSnQtD2aBKRFkDITdJpDFk0KOvuVJgEJ\nEaScRJgtpLxJESaR0yQgFYKUjWiT7tQ2KUKNBE0CkiBIeYg8d1Vhk6LVSNAkID6ClIEkU8BVNSly\njQRNAiIjSNolqZFQ0qQkNRI0CYiJIKmWsEYieZMS1kjQJCAagqRX8hqJhE1KXiNBk4A4CJJSSmok\nkjRJSY0ETQIiIEgaqaqRiNwkVTUSNAkIjSCpo7BGIlqTFNZI0CQgKIKki9oaiQhNUlsjQZOAcAiS\nIsprJII2SXmNBE0CAiFIWmRRIxGoSVnUSNAkIASCpEJGNRLem5RRjQRNArwjSOllVyPhsUnZ1UjQ\nJMAvgpRYpjUSXpqUaY0ETQI8IkgpZV0j0bNJWddI0CTAF4KUjIEaCecmGaiRoEmAFwQpDTM1Eg5N\nMlMjQZOA/ghSAsZqJDo1yViNBE0CeiJIsZmskWjZJJM1EjQJ6IMgRWW4RuJgkwzXSNAkwBlBisd8\njcSeJpmvkaBJgBuCFEkhNRI7m1RIjQRNAhwQpBiKqpH4q0kPT1VhNRI0CejqXeoFFOHux/sPv7xf\nzZrUC4nt7OOppGhU1RfX4yS/gXXzuG6ekrz08GQwGtfTm/nFz3+P/+pAdghSDMOT4+nNPPUqUlo/\nPKb6DawfHquqKvz3D2SBIMVQ1Jd1L91d3Y+qev3wODw5/nD7flgP4q9hMVmuZk2qbwu/fD8Zjesk\nLw1kh2dICEW+rLu4HkuNIsw+1+bL95OL6/Ho/CT1QoA8ECQEsbWLIcLsc23+rBEfj4C2CBL827mn\nrqgmUSPAAUGCZ3t2eBfSJGoEuCFI8OngeSPzTaJGgDOCBG9ann413CRqBPRBkOBHp7sYTDaJGgE9\nESR44HAzkLEmUSOgP4KEvpzvqTPTJGoEeEGQ0EvPW1MNNIkaAb4QJLjzcod31k2iRoBHBAmOPE6U\nyLRJ1AjwiyDBhff5Rtk1iRoB3hEkdBZo2l5GTaJGQAgECd0Enf2aRZOoERAIQUIHESaRK28SNQLC\nIUhoK0KNhNomUSMgKIKEVqLVSChsEjUCQiNIOCxyjYSqJlEjIAKChAOS1EgoaRI1AuIgSNgnYY1E\n8iZRIyAagoQ3Ja+RSNgkagTERJCwm5IaiSRNokZAZAQJO6iqkYjcJGoExEeQsE1hjUS0JlEjIAmC\nhH+jtkYiQpOoEZAKQcJflNdIBG0SNQISIkj4QxY1EoGaRI2AtAgSqiqrGgnvTaJGQHIECfnVSHhs\nEjUCNCBIpcu0RsJLk6gRoARBKlrWNRI9m0SNAD0IUrkM1Eg4N4kaAaoQpEKZqZFwaBI1ArQhSCUy\nViPRqUnUCFCIIBXHZI1EyyZRI0AnglQWwzUSB5tEjQC1CFJBzNdI7GkSNQI0I0ilKKRGYmeTqBGg\nHEEqQlE1Ei+a9FhRIyAHBMm+AmskpEmLyXIxWVIjQL93qReAsIqtkRjWg+HJ8fDkeFgfp14LgAP4\nhGRZ4TWq/vymLubscwDOCJJZ1Gjz3Cja7HMAfRAkm6jR1i4GmgToR5AMokY799TRJEA5gmQNNdqz\nw5smAZoRJFOo0cHzRjQJUIsg2UGNWp5+pUmATgTJCGrU6S4GmgQoRJAsoEYONwPRJEAbgpQ9auR8\nTx1NAlQhSHmjRj1vTaVJgB4EKWPUyMsd3jQJUIIg5YoaeZwoQZMADQhSlqiR9/lGNAlIjiDlhxoF\nmrZHk4C0CFJmqFHQ2a80CUiIIOWEGkWYRE6TgFQIUjaoUYQaCZoEJEGQ8kCNotVI0CQgPoKUAWoU\nuUaCJgGRESTtqFGSGgmaBMREkFSjRglrJGgSEA1B0osaJa+RoElAHARJKWqkpEaCJgERECSNqJGq\nGgmaBIRGkNShRgprJGgSEBRB0oUaqa2RoElAOARJEWqkvEaCJgGBECQtqFEWNRI0CQiBIKlAjTKq\nkaBJgHcEKT1qlF2NBE0C/CJIiVGjTGskaBLgEUFKiRplXSNBkwBfCFIy1MhArdmuMwAAAXpJREFU\njQRNArwgSGlQIzM1EjQJ6I8gJUCNjNVI0CSgJ4IUGzUyWSNBk4A+CFJU1MhwjQRNApwRpHiokfka\nCZoEuCFIkVCjQmokaBLggCDFQI2KqpGgSUBX71IvoAiLyfLiejy9maVeSALr5mk1b0bjejVvVvMm\n9XJiG9aDxWT54ef/Tr0QIANHz8/Pqddg3HQ6/f3331OvAin99NNPqZcAZIAgAQBU4BkSAEAFggQA\nUIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQgSABAFQgSAAAFQgSAEAFggQA\nUIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQgSABAFQgSAAAFQgSAEAFggQA\nUIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAVCBIAQAWCBABQgSABAFQgSAAAFQgSAEAFggQA\nUIEgAQBUIEgAABUIEgBABYIEAFCBIAEAVCBIAAAV/h90BlRfclZCBAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%% Lshape adaptive grids\n",
    "mesh.shape = 'Lshape';\n",
    "mesh.type = 'adaptive';\n",
    "mesh.size = 2e4;\n",
    "pde = 'Poisson';\n",
    "fem = 'P1';\n",
    "% get the matrix\n",
    "[eqn,T] = getfemmatrix(mesh,pde,fem);\n",
    "% compare solvers\n",
    "tic; disp('Direct solver'); x1 = eqn.A\\eqn.b; toc;\n",
    "tic; x2 = mg(eqn.A,eqn.b,T.elem); toc;\n",
    "tic; x3 = amg(eqn.A,eqn.b); toc;\n",
    "fprintf('Difference between direct and mg, amg solvers %0.2g, %0.2g \\n',...\n",
    "         norm(x1-x2)/norm(eqn.b),norm(x1-x3)/norm(eqn.b));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The finest mesh is several uniform refinement of an adaptive mesh shown above. Now the multigrid outperforms the direct solver around the size of 7e5. amg is 4-5 times slower."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: 3-D Linear Element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Direct solver\n",
      "Elapsed time is 0.443978 seconds.\n",
      "\n",
      " Multigrid V-cycle Preconditioner with Conjugate Gradient Method\n",
      "#dof:    35937,  #nnz:   202771, smoothing: (1,1), iter: 11,   err = 5.50e-09,   time = 0.18 s\n",
      "Elapsed time is 0.150138 seconds.\n",
      "\n",
      " Algebraic Multigrid W-cycle Preconditioner with Conjugate Gradient Method\n",
      "  nnz/N: 5.81,   level:  4,   coarse grid 305,   nnz/Nc 30.80\n",
      "#dof:   35937,    iter: 10,   err = 3.5316e-09,   time = 0.92 s\n",
      " \n",
      "Elapsed time is 0.645183 seconds.\n",
      "Difference between direct and mg, amg solvers 1e-09, 4.2e-09 \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4gsWBQUo2/O3LwAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAyMS1Ob3YtMjAxOCAyMTowNTo0MBueoAYAACAA\nSURBVHic7d3Nkd02t4VhtusLoD1zBCoF0WGoetxheOChBwpDY1WHoRRcpVIEKk+sDPoOcE1DJA8O\niL+9Nvg+dQe++qRunD8sLhCHfHh7e1sAALD2i/UAAABYFgIJACCCQAIASCCQAAASCCQAgAQCCQAg\ngUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAAABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAg\ngUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAAABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAg\n4X/WAwBw7OHhYf3vt7c3w5EAY9CQAC0P//r4/Y9lWV5enz9+/yMOJ2BWNCRAwho5IYeWZfn65dv6\nH2smUZUwMQIJMBaSZs2h4OuXb++f3sV/Ev7Cw8MDmYRZEUiAjX0luuX907s1n6hKmBiBBIx2WIli\nIX5+/+3Pj9//WBfuVlQlzIpAAgbJrET7xbrl55K0/hCqEiZDIAF95S/NnUVVwmQIJKCXu0tze2sN\nCut18f+0L0kBVQnTIJCAxoor0WHe5Fir0kIswTMCCWimoBKdcqskBazgwTsCCajV5CxRcT3aYAUP\nfhFIQKGGuxU2abQ/gbRKl6R4PFQluEMgAaf1XpprgqoEdwgkIFenDdxnF+tySlJAVYIvBBJwX79K\n1OrUUQJVCV4QSMBN/b7TmpA4gbTKL0kBVQkuEEjA1rAcGlCPYlQliCOQgP+M3K1QmUZnS1JAVYIy\nAgkwWJob3I02qErQRCDhukxOEaXlnEBalZWkgKoEQQQSrsj2i0S29ShGVYIUAgkXolCJ2qZRTUkK\nuDArdBBIuAQX11YwxAoeFBBImJlCJYql28ypE0ir+pK0YgUPtggkTEgthwKdU0cJVCUYIpAwFZbm\nmqAqwQSBhBloVqLY3XpUtl4XNFy1W1GVMB6BBN9cVCIXi3WHqEoYiUCCS/qVaLAeJSmgKmEYAgme\neMwhv/UoRlXCAAQSfHCxNLeXmUY1J5BW/UpSQFVCbwQSpHmsRHOjKqEfAgminFai2PjFut4lKaAq\noZNfrAcA/OThXy+vz8uy/P7bn7//9qf1oErMceooIVSltcIC9WhIkLBZmvv65dunD59fXp/DnB4y\nyXVbSmhyAmk1piQFVCW0RSDBWM7SXPhf207cXU1fj2KcVUIrBBJsJHYr3JrNP37/w0VVMk+jkSUp\n4B4WaIJAwmjpSpSeSd1VpUthBQ+V2NSAQdbdCh+//3E3TtKRE6qS5maHs9WkU7iGktT8x+ZgswOK\n0ZDQ16kvEuXP5ppVyXyxTgRVCWUIJPRy9otEBbO5l7NKJsafSdpgswPOIpDQWP21FU71Hp2qRD3a\noyrhFAIJzdRcW6FyNjevSmXj752j5iVpRVVCDgIJteorUZNJU6cqYbUeJfz+258vr89UJaQRSCik\nedlTk6ok0kIOmZSkdQPk5lXgrBLSCCSc1vayp5vpsr7iDK5KxdP9fGXuVg4FIRo5q4QEAgm5elSi\nfgfv5meVdAwoSemnOqzXxX9CVcIhAgl3aC7N5RhQlZQX6wZIV6K9OBqpStgjkHBT7zsSjZnN+1Ul\nR2nUvCS1ekqpSogRSNgaU4kO58dObWatSku3U19nOT2BdLYSLf8+0vU6RvtopCphRSDhPxPcpDXh\n4vvCa0pSQQ6dRVXCQiBhsThLZLjY1WQFz9FiXaUeq523opF7WIBAui6r3Qrms3llVTIf/wADKtEt\nrOBdGYF0RZpLc4MX0wz3hVstG+as2rV9Tm490rsjYQXvmgikC1HYwC1VLwqqktT4GzKsRLdQlS6I\nQLoEkUqkOZvnVyXN8Z+yqSaGOZS5yYKqdCkE0swUKtFKeTa/4AY8R5exoCpdB4E0IakcyqQQBumq\n1CpQbR/pyEp095Ge2olOVboCAmkqIktze8r1KHarKnkZf0KctevXVH2hKk2PQJqBeCVyN5vPdGHW\nw0qkc+O+gpFQlSZGIPkmW4lOUViv24irksjcfYrgrrmGqEqzIpBcEq9EMY+z+ap5VRoQvZkD7l2S\n8h/p+6d3P/758fjr49lfQVWaD4HkiaMcClynUfDy+vz+6Z1gh9twXYm+//X349PpQFqoStMhkHyY\nY2nOnfjmPbJnlYoHpnMmqRJVaRoEkjR3lSiWOdnJlo/DuyTUxFLbR+q6Em3URyNVaQ4EkijvlWiO\nQ+8986/Qts2hTiXJ6vmhKnlHIEmTbQ/Tu3vpz/EreLJrhk20ikbuYeEagSQtTHzu5iDv9SjzGmvL\nkCOGTx8+x7+xuWnOJMVYwXOKQFKnfDr90KnZzWPcxvJfnYJHGn5y2OZXPMLxil/T5tHICp47BJKu\ntR6Zn7e4lIILByxNX53NWaIBl/mZsiQFVCVfCCQ3XFQl7/Na8fjrX52Zds0V6xSNVCUvCCRpm3NI\n4lXJexpVKn51FI4zJi5JAVXJBQLJHxdVKYdasrba5bV/dQ4f6ZSVqP417RqNVCVxBJK6w412glXJ\n+/F1w/GnXx3ZHJq+JAVUJWUEkmM6VekKE9lZ+1dH5MXSNyAaqUqaCCQHEt9GEqxKHnWa/uKrDbl4\njS5SkgKqkiACaQa2ValsCtOZoPtNwevq3EU0fE2HRSNVSQqB5MPdSzbUX/qzjPcD6k5Xcgv/semv\nOhmccKmSFFCVdBBIU3E08c3n7m4FnXN+LgyORqqSAgLJjfzr2g2b+IrnC5HIbDXf5T/bVkX2lAuW\npIALs5ojkOY0oCp5n7Pqx5+zgfvwJZiyyPZ4OCbRyAqeIQLJk7MX/2aNqJMmzyqvjjJW8EwQSJPr\ndDB+zXrU/DutslVJZ9XOcCRUpfEIJGfK7pDU9mC8coIwn3/Pjr84h9TO+aEAVWkkAukqZA/GlQ3L\nCcFX51Q16Tpy87pGVRqGQPKnrCTF/3apmGRFVnKK5Yzf6nJzVCVlVKUBCKTLqTkYnz6NzPNAqiqZ\nVxOdkVztohtWCCSXakpS/BOWM5NvkxlBZKrd6FGJah6pr6qk+Zq2Er8QxFJvBNJ1SR2MD7APVNk7\nQSzRq/Py+mw4DPNqEo/kxz8/Hn99HPPr9u+N8Ekhk7oikLyqL0nxz1nuTcois1Kxzfi99A9fVam3\n73/9/fjUPZAOn/DrHLfZIpBwvyp5T6OVciW65eX1+f3TO8MJ8SIlyeN7Yz4EkmMfv//R9j6nvQ/G\nrWbVr1++ffrweRk41zR/pLJVafBr2rwkFV//CT0QSL59+vC54Ufl8NKfIkfHZWY67DW8MKtOSWoo\n85kkjUYikLAVr+A5nYbiHHL6EG652laUjfponOkYZT4Ekm+ttjbc+smLt8/tNPXu65dv8X9vHsX4\nV2eCklTwjF02+K0QSLip7en0rp/tw8Neqwm05pGuORRGvv6/4T/ih6NQlUx++9lo7H0pQjREILnX\nqSStn3nxqqQ8tlP2kRNegvjPbauSr5LE0pxHBBLuUzgY37g73XiZOjeVKP7z/fg3+RQIvjq9paOx\nSUJf6vnUQSDNoHlJOvy01xyMN1z3WweT+Gsu0mgfLbfE869tVVIuSQ0rEWlkhUDCVmLGMTwY97I0\nd/fJuVWJNn8nPekrVCXbWXuNRpbmZkIgTaLfdrtbv2sZMgUUTDeah/A5ObT+zVuLdfGfW1UlnZLU\n6cvO1CNDBBJ+kjnXDDgYL5tVRebKWP7S3HJ+/ApVabDNl8ya//ApnzQvCKR51Jeks7Nh5sH4qVGF\nH/jy+txjuuktfqT5lSjTrWoyviqtKTiYl2VbFCOQUKXVwXiTMwEi9ehUJdr8w+LxD65Kba9ZlZZ4\nb7RdP6QemSOQplJTkmo+2JUb8Mr+4YZ5GoU8KB5Gzj9Mz7+Dq1LvJ3zwbgXSSAGBhGVpMbmcvfTn\nNJujNktz5hVtgrNKp+JTZ5MF6hFIsxm53e7wty8/T3ybwXTKIZMpaTPvV14xKHP8OfNv76q0PtK2\nT7vhMYqXqJ4egYT2s/nhxNfvjPTgNGq+W6HT+NNVaVHqppXjqSxJpJEOAmlCp0pSp9kwnvi6Tn8j\n06h4t0Jb+fNvfIVWwRW8aZZt0QqBhI7EL8yaqXkl2vzw3gnXbwWvrJr0yKHikmSeyogRSHPKLElj\n6kW/01pdx38qh8oeYPH4C+bfVpsdal5KtaMT0kgNgXRdwxa7wse++WTUb/wiS3PN9ahKOdE4ZmmO\n7XYTIJCmZbvdbk/kvEVC16W5w19X84uK59+R+8LXnBOMCuW34mURSBc1coKIo7FVVWo4/sE5tP5S\nwwm6bVXaR6PVboX8kCaNNBFIM7tVkmxnw/qD8Vbjb7U0ZzW7VS5SFVSl9CNl1xwqEUiXY5JG+2g0\n3IBnUok2AxBZv2pVlTrdCaJA5jkthaFij0CanNqZpFhZVaqZzRV2K7RNoyZn8ovPKnW9E0Qnsh8H\nLATS1Rgem9+KxlMH42Xj71qJ5pjgzlal/R/qbHJ7//Tuxz8/Hn993P9PX798m+DFmhiBNL81CUTm\ni71OW7zMl+b2erwEDZPg7tWGFr3vEh36/tffj08HgTTyrhkoQCBhnPT64d2qlD/tKizN7ckeEMQO\nq1LmbgWdknRoji47NwLpEpTrUSxx6c+c8QtWojGaJ8FaldbdCmsmuXC4E5000kcgXcLXL99EFity\nNlmssfTy+pz5Y60qkflFbDtZ4yf/JQjESxLEEUgQFa/g3ZrjvFSiMVdQbfJbbm1hcF2SqEdeEEjz\nCx9Lnf3f+SN5eX1+//Ru/5e95JAjibNE646YzB+lVpJE3vbIQSBNTmpqKPPy+hxXpcVbDg17CQqS\n4O5uhXWf9K194crWM2HwgkC6EI8lKcyA4UyGyOBXOjf4KJC5e3tz6vFwX/ieTkkSOXWKTATSzEQm\nhUrro5jjdn+2Tl1u7jBxHVUl8a/fYY9Amtbh59BLSbp1lkj/Hhax8VNhopqczfL0RQ3uViWdkgRH\nCCQIiXMosfrvoiqJTMfFV+C+u9glXpXWoxai0RECaU6JT6BUSVrHuTncDn+ezqTFtCqlf7XhDBie\nt/CF1qU0s/Of2ERVMkwCkXc4ziKQJuToePDTh89hw0LZgF1UpcHqn5Czs7l4VVooSX4QSFekUJLC\nYfXhOedTc4d5Vdozmfs2S3NWp6+Wn48tTJJA6s2AUwik2YgfCe53K2zOVZSNX6cqDX7+m9+ktWY2\nV6hKt8ZPSXKBQJpK/kducEkacG2FkVVJ4Rg8HcBl82+Tx7WpSnw7FfkIJHSXnhnjaKw/hrWtSgOO\nwZtXoh6sqlI6UClJ+gikeZz9sHUtSfFBcf5NjJpMFol7WHTVe7I7+4jOzr/N3wxxVRpQkhRqKyoR\nSJPQOfQruNxciMazdzq4+zOXPpPU4ImvcgN3pk4PKt7Tr4CSJI5AurSGJan+LFGPaWLYCl7zaS6+\nI1HNU6ow/4ZhdB0J9WgOBNIMDCedJrsVwlVqOs0pAzY7tH3+Nwk602JXp4u1nxq/SEjjEIHkXuWn\nqzgJWk0uY2YHnX3ht3TarXB3/h3ZLRT2hUMZgYRz+m3g7r0TvUlVOrxbYOVTIZ6UrazRmHkPi0wF\nLyglSRaB5FuTz1VmEvRYbxk/L7StSjXjH7aBO30JcJMgbFWVOHU0GQLJsTGzeb9KtB//mK/r2l5t\nSOeLROOfgU00tq1KNSOBCAIJy3KUBAOurWCrviqV3S/cqpFozr81VYl6NB8Cyat+88uYg9Zb4x95\nTaOzVSn+m/nPv04lihku1u2fuoKqVD/+90/vfvzz4/HXx5ofgrYIJJd6pNHL63P4duqYZUCdo/V+\nG/CkdivESSDYLU5VpfTdbPN9/+vvxycCSQiB5E/b2XyzNKeQE+PvjnG2Kt3dSB3/WGwk1g8zq9Ld\nu9nmj4SSJIVAuq79J//w7kQ9fq9C7O1lVqXEjrX15/QYXr33/95MVnaES0ZVanuwQkmSQiA5Uz+b\np3cr9J6tMsff4+p2ORIXZk3Mg1JLc2kiaXR3k8WtqiS42IiGCCRPKtPIaoutO4kVvPgl0K9Eh2Qb\n6saYyzrI7j+8JgJpfqc2cPe+J0X+J9/8Puv7Fbx4U8DiLYeWfyNW5NrbmUkQVyXq0fQIJDfOHsep\nfZHI43HopiqNuRNEJ+ujcNcJ1oOATku47p6QiRFIPpz6wFQuzZlXk1XYiW4+Ev0Lszp1NglIjukR\nSPNQq0SxCeYR12mkkOs1NvVuaf0mJ+pE/GI9ANx396Oy3v0s/F/9b1w7QRM1H/W2I7mmfRq9H3JP\n8Rw5I9mMP75eOCZDQ1KXmM2VK9FMXM99rS5qoKZ5VaIkKaAh+RP6UNtKtNewmlSO0HZjmPdJat2I\nsaFTMtIjSSw2UpXmQ0OStpkN3X2RqNVkIfJ1TnemOXWU0LAqUZLM0ZB0rZ+NTSUaNgCd8zdWJcn1\n9HR3sU6nXlSOhKo0DQJJ2oCluX7azua3lp76cZ1Gi8Uz1tbZehcyqTKWCDZbBJI08xwqLkneZ3Pv\n48+czXXm381IyhYbqUreEUi6Pn34LLJipkBn/VCf91NHlSqrEpFmiECSFmZh24m4IAm81wvv4z9F\nZ/5dR9LkbrA6jwv5CCRd6yV8fJWDfrP5mOfBexp5r0cNx19clQgzKwSSD7ZVyVciXlnxqZdZ59+1\nKs36ACdDIEmLk8BFVepdL3o/A97rkXefPnzucUnvghW8iUNaGYHkjFVVykkC77O59/HXLHYpzL+9\nFxupSvoIJH80q9Kw2VzwsSvwfupo1TUaT1UlhZC+GgJJ3a35d3xVmjsJvNejerbz78hApSrJIpAc\n06lKg2fz5o/aexp5r0f7G0z0TovMqkRJGoxAciA9/46sSocj8T6be9d8n3STH+UCVUkNgTQDnao0\nTMPHS6DaOgzUYdF4typdLaRtEUg+5My/Y6rSZiTeZ3Pv42++WDd4/hVZbKQqiSCQpjK4KtnO5lcr\nhXsis3kng6MxUZXeP7378c+PYSO5MgLJjfz5t3dVsr2Fa0Pe61Enw5JAMFBvVaXvf/1tMp6rIZDm\n1LsqffrwWWE2r3mMCuOvITibn3LqbrAjcWFWQwSSJ2fn335Vqcf1XZBvzEUN+v38u3ezNbepSkTU\nGATS5Naq1DyWRG5IWlaSvNcj7/LfPIZJQFUaj0Bypmz+bbuC53029z7+MYt1/eZiX4uNa1UinAYg\nkC6kSVVaZ3OdTW46IxnA12y+VzB+8ySgKg1DIPlTM/9e8Cu0G97r0UjMwrH3T+9ElqknRiBdUXFV\n2szmOtmWORLvaXTBehSYR6Ph7TEv5X/WA0CJ9e7mNT9hOTlBeJ/NvbtsGhmKb48Z/7/ohIZ0afVn\nlRyVJNeBarVP2ryarMZf0ygkaPi/xWegukND8qq+JK0/Z8n4sLmezRf/4//04bPr2dDL+DeVKP5z\nF+P3jkDCskT14vBTl57NW0VjvXBNI9fBc8j26Q3VpOZZ/f23P5t8k7p+JLfcyiEMRiA51jYJ5liX\nODwS955Srl8RcYnjsPjv8BKMwTkk/GR/VilnNpc6k7Q50+A9jRS2Gtecv2k7m7c6kxTe5OtZovTf\nJI2GoSH51mO5LK5KHmfzuCR5HH+s1WKXFcHZPKcSwQqBhGPps0q3/j6f8yZCCfj04fPL63O/Eyen\nmA9gVfaElJ0l4i09GIHkXr8kCLOhuyPK9QlRmMfPClG0Gfa6TuXr4SjM5jW7FRTGfzUEEo7F16xb\n8j6cUiXJVxrFtzkI/7Fmavzn5g8qfwBd3wnhFq6Pvz6mB7C4OpDCQiDNYUASnF3BsyUVjWmH1efW\n4J1WpR6+//X349NBILXawO3l/TMZAgkHDg+Ec6qSSBKIXFwgYV+JDv/O/n9dq1L633aSef7G5D3Q\n8IBJ4T18TQQStnK+BrsIV6UwfpFo3MjJocx7ey8CK3iHRt6uaT3HuQi/IZGPQJrEyPk3XZU0k8Bc\nZqc59dSZrOCJ7Prrl0O8ew0RSPhJ/lyjWZXi8StEY04lqqFWlQY84fG7rvnarPkb5uIIpHnUz79n\n57VbVckqCXTm5aXoNE/8pJ16LIOr0q2S1PVFH7A0Z3VJdawIJPy/4tlcsyoF46OxuBJVjlOtKrWS\nzqG264deLkk+MQJpKlbVZF+VTJLAcC7uvTSXaVhV2idB85d78FEOi3UKCCQsS6PZ3LAqJcbfOxqb\nBEDDEZpUpYbjP7s016QkkUYiCKTZFMy/DSevuCrpXAK8h4aVaP961b8iA6pS8+12sgu/GIZAQnuD\no+jutNi2JHm5VsKwqlT53NbvVqiMRuqRDgJpQqfm304T1lqVmv/kjWFrU53OEvWeDbtWpfC91LLx\ni3yhlTSSQiBdWu/ZfL3dn/lnvrgkdd2tMOyiBkuf1/rrl28Ft2vqsTQn8nVdVCKQ5qTwndBVv80O\nXecgq6W5Tg+qR1UK+6QzByxSiWI6nxEEBNJ1DTuiXD/2bT//BV/jzRnAsA3c42fDthdmzR//mN0K\nZ0sSaSSIQJpWev41Wd9Q/grtMrYSGc6GTVbw4vEnLtwQ/iO+sROQQCChuzgaW1Wlsvn0MKRFvtMa\nDDtQ6LfZwXBpLr8kUY80EUgzu1WSzE//VlalJuM3zCGR2bC4Kt0av3L9jYk8/9gjkC7HcLEungV6\nnFXKH8nL67NVJKvNhmer0n78UrsV2G7nGoE0OantdnsFVal4uolPYzBhxWqq0tLtThCdKH8cQCBd\ni+HB461oPFWVCsa/X5rL36ncXPph2h7a51SlMP7DSqRTTRIjIY3EEUjzW5NAZL441GMDXmJ6NbnR\ngP5smFOV9B8F/CKQME56/fBuVcoM1JzdCuPj2dHN3/ZVKfMskXhJIkr1EUiXIH4mKXarKuXMdKdO\nzg8uSb5u/haewziHvF+43cv7/+IIpEvQOeGcE41rVcq8TprUF4kOZV4hQmf86zFBeG43p44SxEsS\nxBFIV+GoJAVxVTqcWSpzaNgT4mixbr80N8eE7uudf2UE0vzUjhPzkyB8W+jWtRWkHlSC/mJd+hTR\nuiMm86fpVJMwEv3nHysCaXLx1OCuJC3/1oswY64reA3vb9v7CRF/wu/ubFzHb3JndFwNgQQDZ5Mg\nRNGnD58XjcsBZDK/TeItxddWyLysg05Joh75QiDNbD8pOC1J4VH0uNqQxyekxqkvex0+M23vYdGV\n/tfvsEEgTUv8c5hIgs3pisNMXeSrklTOFVSi9PjvruDplCQ4QiBBSHzcvU5n+4PxtlWpR0kSSaPe\nlz3tdw+LevHZL6LRCwJpTolPoM4iVXxNo/Anh2O+dTDupSrd1WO6rHxm8t8hiapkmAQi73CcRSBN\nyNfxYOZoDw/GW1WltiFtOBs2qUQF41euSgslyQ8CaTY5HzzzkrRWolMjcVGVrJ5Y82fg8NUxSQLq\nkV8EEsbJv7bC3bPlS+uqZBXSlfN187NElU+CeVW6NX5KkgsE0lTyP3KD59/EJFUwEtmqNOwp7bRb\nocn4N6/OGlEDOLpKEw4RSPMQPAAMM1GngTWvSpUhPSaNzJfmMplUpfTXYClJ+gik6+pakuLJKOec\nVtlMIVuVmuu9gXvpEKjxqzOgJHHqaAIE0iREDv2Kr8AdH9uefSwNq1JxSJf9q7uPdEAOrb+o93eV\nusocPyVJHIE0g+LPWKuSVH9HospLvNhWpR6z+TQNb31vkAS4i0BClYbnCeqvg5muSkve/G67J35Y\nJdr80h6/bn/djU5nlU6Nn5KkjEByr/LTVTb/yt6kdczVhmKtfqZVJWr+nCTeGz3uYcGpo5kQSL6N\nP9brt3UqRGO4KV/lj6pcwcsP6ZrZMAzPpBLFY2j1exM5dPht2cO/OQYlSRaBhKz5V7YSJXS92lAl\n2xxaNbldUEG6tLqHhfnriLYIJMcGHOUNzqHm52+Kq1LOSIr34y3VmziaqHyq698blSt4NeOnJGki\nkLxq+3Haz78miyqd9gf3qEpn/6FIJVoVP/DmxyjmK3jQQSDhJ+ZLc502uRVUpVYjOfz5Tg/PC4ad\n808KqlL9S0NJEkQgudTjg/Ty+jz9inyrqpTzl9UqUezUg41ra9fpO78qTf9GvSwCyZ+2aWReifa6\nfhPoVFU6PNOTHptyDgX5z+34lbQe+8LTv46SJIVAuqjDHLL9TujIqSG/Km32oSX2Sed/kUh8EjQ/\nRklXJerRxAgkZ+rnMhcnkAdEY2ZV+hjdYX052ietX4liiWe1Rw4Vv11vvTo9rgArfnxwKQSSJzWf\nnMzpxrYkjZdTlT59+Bz+fPPMuLvc3K1XVvYYhQ14V0MgTc58+aXYsGi8W5WCdbHOVyVKcPHeiF+d\nTu8HSpIOAsmNs5+Z4kNLk5JkPiPcvTDr2pMqnxmrRxq/pmNqR8NHGtLo5fW5yU+DLALJh/zPtovD\n3kyDozFdlVz3oTB+7++NflXm/dO7H//8ePz1sflPxikE0jwaHvZOMAUXm++8RdcbyY+xvhv7vTrf\n//r78YlAMkYgOZCeSrwf9qaZRGOrS3+aCw+hyRVUDW3eANO8OtgjkNTdSqOZckjz4L3HlzTHPNL4\nvWHVdHs/0uavDlsbFBBI/ow5NtRZtbMdybpG5MJ8vSH90s+3vnpxBJK0+JBtpkrky+Dr2RS49d4Q\nOaQoljP+hq8OJckcgaRr/WxYHQNSkmKCB+PpYxTzZ2wkwVcHBQgkaVf4jDk6Jq08GG/4SMXfGE0e\n6dlAbVKVKEm2frEeAG769OHz+6d3tp+NzdUKDOmMJMxZJieWwu8NM2b6jeG9HhWP3/DVQT0akqi3\nt7dlWR4eHlxPK7Maf1bpVCXynkaJS6rnqHx1KEmGaEjS3t7efv/tT9tm0LWanPrk65SkYMDBeH4l\nmsl6laYaVCWPaEjqqErK8g/GT0VvzY5K83pUWS8ajr+4KlGSrBBIPry9vT08PCxGV5hW2OS2jkRw\npmi4xavy54i8TMV6jJ8NeI4QSG5QlQLNC+HUb8CLfw7aKnh1KEkmCCRnrKpStq7gEwAAB3lJREFU\nj5I03wf+1sF44pE2PHinHqVRlfQRSP6sVWnxf4+4Ajrrh4cyv8vcvBKJPCfFBxm2N2O89ZfnO2YS\nRyB5NX4FTzwJpNya+FiaU0BVkkUg+Wa72cGKl2iML8zadQZ08WwkyN5hhJI0GN9Dcu/t7S18XWnA\n72r4TaCrfc77fZFIJ43KXlPbS7n7upr79AikSSh8hXYktS/J7q1faF16fkmz8qIGWO69OiTWSATS\nPMZUJf0ksHXr2gqdDsabXNTAkEi9oyqJ4BzSbK5zVknqTFLmboW2p9N1Hn4ZtfHfenU4kzQMgTSh\n3hvw6pNgpo93OmD2j7TVhVnVZvM5XlP9mzHOjUCa1hWqkmFJqtzAzc5jtUCN7V8dStIYBNLM+lUl\nqeWykRp+kajmYNz7k68/fqqSCQJpfnNXpWHR2KnQFFQl/dl8GvGrw66HAQikS+hRlYqTwNdRZ00l\nyr9ETf5f1nRq8L4CdYJXxxEC6UJmrUqdStLgczyZVcnXbL7ndPzUozEIpGtpe2HWKc8kGV5u7u7B\n+HzPtiPvn97xDbzeCKQrMry1UqeljyY70cN/mK/MTLwBz2mgrjkUPjjoh0C6riYreJq3cD2l0+xf\n87QcXvpTdjbPfKSy408IUUQODUMgXVqTqiRyC9ezJUmnEt0Sr+B5nM39ohJZIZAw7WaHW3wtiM1x\nOt1FoJJD5ggkLEtdVcqvJr0X99IjGVmJ2j7S0EH9Lo3qpxFLcyIIJPxnyqqkvzSXts7mmpsd/Mbk\nQiXSQyDhJ2VVSWf/dzwSwem7kscvaYq8MTaoRJoIJBzwXpWmiaLD2VyzKh1SSyMqkTgCCcfOVqW7\nJWnAcX2Yps3rWqtHmngUIlXJfACZyCEvCCSkuKhK3s8SFROvSgr1iKU5Xwgk3JFflQZXk0QOud6T\nFmQ+kyJVac82jahEThFIyCJVlXJqgcjXdcucnc3VqtLXL9+snnwqkWsEEnLlVKVbJanJIfyppTmr\nkmRVVsZXpcTvGn80QCWaA4GEc0yqUtnhv9OSVLPYpVCVRi7WkUOTIZBwWvoeFg3PJF1wt0L9U2d7\nVmlYGrE0NyUCCYX63cOiVQ6Z7/82pFCVeqASzY1AQpXDFbxNEuQfrXufQyt7Sdv4PLyHRSuHj7Rr\n/FOJroBAQq36qtRvac5RSeo0zmEreJ3GTyW6FAIJbWyqUmYSeK9ErfTeJ+1uBY8cuiYCCc3cqkr7\nSXDkbgUXJWnAhsCuVanhM8zS3JURSGgsrkrrcW4w96654rl+ZF42qUqbR9pk/FQiLAQSeoj3hQe2\n60XKJWn8wKSuNkQOIUYgoZe1KonMfYi1OqtUHKgszWGPQEJH/b6rdFbvklQWura9rb4qFYyfSoQE\nAgndSV2YVYfIKuLZqlRzqmwhh5BEIGEEhaqkfCbJVllVynwyqUTIRyBhHKrSSjAaT1Wlu+Mnh1CA\nQMJQtlWpU0nq1C3Gu1uVch4pS3MoRiDBAFVJ2d2qdBioVCLUI5BgI30Pi37MzyTJ1qNYoirtx08l\nQisEEiwpbHYYyUUardJViUqE5ggk2Bu8gte2JM39td/4HhbLsnz68HmhEqEbAgkSrlCVfNWjWJy4\n5BD6IZAgZFhVGn8myW8asTSHYQgkaLlCVfKCpTkMRiBB0YCq1KQkZZ5A8lWPqESwQiBB1DRVyUsa\nkUMwRyBBWteqZP6dJBEszUEEgQR1rquScuBRiaCGQIIPnarSx+9/FH+R6O4/lE0jKhE0EUhwo1NV\n+vThs2ZsNEclgjgCCc54uTCrTj0ih+AFgQR/2lalsq0N6fU6kTRiaQ6+EEjwyktVGo9KBKcIJDjW\n6h4Wbfd/G9YjKhFcI5DgntS+cJM0ohJhDgQSJlG5gneqJInccoIcwmQIJMzDvCoNq0cszWFKBBJm\nU1yVKs8kDUgjKhHmRiBhQuOr0tcv3/r9LnIIF0EgYVoFVSmnJB2eQOp0uQeW5nApBBJmNqYqNV+s\noxLhmggkzO9UVTp7JqltGlGJcGUEEi7BfANeGpUIWAgkXEpmVUqUpM0JpMp6RA4BMQIJ19KwKlXu\nEV/IIeBnBBKu6G5V6nR3cyoRkEAg4aIqL8x6Nq6oRMBdBBIuLbGCty9J6wmk/DSiEgH5CCSg/a2V\nyCGgAIEELMuNqnR4Jildj1iaA4oRSMB/ElUprNfdSiMqEVCPQAJ+sqlKoSQl/j6VCGiFQAIOHFal\nuB5RiYDmHvg4AQmhKoX4idsSHxygOQIJuCNUpYDPC9APgQQAkPCL9QAAAFgWAgkAIIJAAgBIIJAA\nABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAggUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAA\nABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAggUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAA\nABIIJACABAIJACCBQAIASCCQAAASCCQAgAQCCQAggUACAEggkAAAEggkAIAEAgkAIIFAAgBIIJAA\nABIIJACABAIJACCBQAIASCCQAAASCCQAgIT/A/atmRzsFwPsAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "% Cube uniform grids\n",
    "mesh.shape = 'cube';\n",
    "mesh.type = 'uniform';\n",
    "mesh.size = 2e4;\n",
    "pde = 'Poisson';\n",
    "fem = 'P1';\n",
    "% get the matrix\n",
    "[eqn,T] = getfemmatrix3(mesh,pde,fem);\n",
    "% compare solvers\n",
    "tic; disp('Direct solver'); x1 = eqn.A\\eqn.b; toc;\n",
    "tic; x2 = mg(eqn.A,eqn.b,T.elem); toc;\n",
    "tic; x3 = amg(eqn.A,eqn.b); toc;\n",
    "fprintf('Difference between direct and mg, amg solvers %0.2g, %0.2g \\n',...\n",
    "         norm(x1-x2)/norm(eqn.b),norm(x1-x3)/norm(eqn.b));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For 3-D linear element, `mg` wins at an even smaller size $3.6\\times 10^4$. Again `amg` is 3-4 times slower than `mg`. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-matlab",
   "name": "matlab",
   "version": "0.14.3"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "30px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
